function [channel_features, total_features] = EEG_frequency_feature_extraction(random_samples, fs, nfft, N1, N2)
    % random_samples: Cell array where each cell contains EEG data matrix
    % fs: Sampling rate
    % nfft: Number of FFT points
    % N1: Start frequency of interest
    % N2: End frequency of interest

    cellData = random_samples;
    numChannels = 18;  % Number of channels in your data
    bandWidth = 1; % 频带宽度
    % 定义频带
    freqBands = N1:bandWidth:N2; 
    numBands = numel(freqBands);

    % 初始化单元格数组以存储 FFT 结果
    channel_features = cell(numel(cellData), 1);
    total_features = cell(numel(cellData), 1);

    % 创建频率轴
    frequencyAxis = (0:nfft-1) * (fs / nfft);

    % 循环遍历每个数初始化矩阵以存储当前段的 FFT 结果据段
    for i = 1:numel(cellData)
        % 获取当前数据
        currentData = cellData{i};

        % 只选取18通道
        currentData = currentData(1:numChannels, :);

        % 初始化矩阵以存储当前段的 FFT 结果
        currentBandMatrix = zeros(numChannels, numBands);

        % 循环通道
        for channel = 1:numChannels
            % 获取当前通道
            signal = currentData(channel, :);

            % 计算FFT
            y = fft(signal, nfft);
            freqData = abs(y / nfft);

            % 循环遍历频带
            for j = 1:numBands
                % 当前频带的起始频率
                bandStart = freqBands(j);
                % 当前频带的结束频率（带宽bandWidth）
                bandEnd = bandStart + bandWidth;
                % 创建一个逻辑掩码，用于选择频率轴上落在当前频带范围内的频率
                bandMask = frequencyAxis >= bandStart & frequencyAxis <= bandEnd;
                % 计算当前频带的功率
                % 将频率数据中符合频带掩码的部分提取出来，并计算其平方和
                bandPower = sum(freqData(bandMask).^2);
                currentBandMatrix(channel, j) = bandPower;
            end
        end

        % Store FFT results for current segment in cell arrays
        channel_features{i} = {currentBandMatrix};
        total_features{i} = [];  % Placeholder for additional features if needed
    end
end
%{
自定义FFT:取单边结果
function [P1] = computeSingleSidedFFT(signal, nfft)
   % 确定输入信号的维度7
    [numChannels, ~] = size(signal);
    
    % 初始化输出矩阵
    P1 = zeros(numChannels, nfft/2+1);
    
    % 对每个通道分别进行 FFT 计算
    for channel = 1:numChannels
        % 获取当前通道的信号
        currentSignal = signal(channel, :);
        
        % 计算 FFT
        y = fft(currentSignal, nfft);
        
        % 计算双边谱
        P2 = abs(y / nfft);
        
        % 计算单边谱
        P1(channel, :) = P2(1:nfft/2+1);
        P1(channel, 2:end-1) = 2 * P1(channel, 2:end-1);
    end
end
%}